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Abstract Thermo-elastic behavior of perfect single crystal is considered. The crystal 
is represented as a set of interacting particles (atoms). The approach for determina- 
tion of equivalent continuum values for the discrete system is proposed. Averaging 
of equations of particles' motion and long wave approximation are used in order 
to make link between the discrete system and equivalent continuum. Basic balance 
equations for equivalent continuum are derived from microscopic equations. Macro- 
scopic values such as Piola and Cauchy stress tensors and heat flux are represented 
via microscopic parameters. Connection between the heat flux and temperature is 
discussed. Equation of state in Mie-Gruneisen form connecting Cauchy stress ten- 
sor with deformation gradient and thermal energy is obtained from microscopic 
considerations. 



1 Introduction 

Determination of the relation between parameters of discrete and continuum sys- 
tems is one of the challenging problems for modern physics. There were intensive 
investigations in this area for the last several decades. However the problem is far 
from its final solution. At the beginning the problem was only of a fundamental in- 
terest. However, practical interest is increasing now. The increase is caused by fast 
development of discrete [1] and discrete-continuum [17, 16] methods for simulation 
of mechanical behavior of bodies under mechanical and thermal loadings. Various 
methods for transition from discrete system to equivalent continuum are considered 
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in literature. Long wave approximation is used in [2]. The concept of quasicon- 
tinuum is proposed in [10]. Localization functions are used in [4, 5, 18]. These ap- 
proaches give the opportunity to spread mechanical parameters determined in lattice 
nodes on all volume of the body. Decomposition of motions on slow macroscopic 
and fast thermal ones is used for description of thermal properties. There are dif- 
ferent approaches for decomposition. In papers [4, 5, 18, 22] the decomposition of 
particle velocities is carried out using localization functions. As a result, dependen- 
cies of stress tensor and heat flux on parameters of the discrete system were obtained 
and analyzed. Another approach was proposed in [21]. Fourier transformation was 
used for decomposition of displacements and particle velocities. Different methods 
for decomposition were discussed. It was noted that the result of the decomposi- 
tion is not unique. It should depend on characteristic time and spatial scales of the 
problem. 

The approach based on averaging of equations of motions and application of 
the long wave approximation [2] was proposed in papers [7, 8]. The derivation of 
expressions for stress tensors for ideal crystals was carried out in book [7]. Only pair 
potentials were considered. Thermal motion was neglected. The influence of thermal 
oscillations on mechanical properties was considered in [7, 8] for one-dimensional 
case. The proposed approach gives an opportunity to carry out analytical derivations. 
In particular, the equation of state in Mie-Gruneisen form was obtained in papers [7, 
8, 9]. 

In the present paper a generalization of approaches proposed in [7, 8] for two 
and three-dimensional cases is carried out. The expressions connecting Piola and 
Cauchy stress tensors and heat flux with parameters of the discrete system are de- 
rived. The approach for derivation of constitutive relations for Cauchy stress tensor 
and heat flux is discussed. 



2 Hypotheses 

Let us consider discrete system consisted of particles, which form the infinite ideal 
crystal lattice in cZ-dimensional space (d = 1,2 or 3). Crystals with simple structure 
are investigated only (i.e. crystals that are invariant to translation on any vector 
connecting lattice nodes). For the sake of simplicity let the particles interact via 
pairwise potential of Lennard- Jones type. Generalization of the approach discussed 
in the present paper for the case of multibody interatomic potentials is considered 
in papers [11, 12]. 

Two main principles are used for transition from discrete system to equivalent 
continuum: the long wave approximation [2] and decomposition of particles' mo- 
tions into slow continuum and fast thermal one [8, 21]. First let us focus on de- 
composition. In literature it is carried out using different types of averaging such 
as spatial averaging, time averaging, averaging over phase space or over frequency 
spectrum, etc. It was noted in paper [21] that unique decomposition is impossible 
because rules for a choice of averaging parameters like averaging time, represen- 
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tative volume, etc. do not exist. The only possible rule for these parameters is that 
they should depend on time and spatial scales of the problem being solved. Let us 

denote average and oscillating (thermal) components of physical value / as ^/\ 
and / respectively. Obviously, 

f=(f)+f, f = f-(f)- (1) 



Different expressions for the averaging operator (^J are proposed in literature. The 
following operator was used in paper [8] for one dimensional case 

v 1 rt+T/2 n+A/2 

where is magnitude of physical value / for particle number k. Parameters T 
and A should satisfy the following relations 1 <C A <iV, 7J,,,,, <C T -C T max , where 
T m i„ and T max are minimal and maximal periods of oscillations in the system, N 
is the total number of particles. Obviously these limitations are too weak. Direct 
and inverse Fourier transformations were used for decomposition in paper [21]. The 
direct transformation gives 



fe m dt, (3) 



F(y) = \ °' v<v "»" (4) 

0, V > V cu toff [ F(V), V > V cu toff 

Here F is Fourier transform of value /; i is imaginary unit; Vcutoff is cut-off fre- 
quency, which should be taken in the range 0.5-50 THz [2]. Inverse Fourier trans- 
formation was used in order to obtain (/) and / 




In the framework of the given approach the choice of cut-off frequency is almost 
arbitrary. In papers [4, 5, 22] the following relations were used for decomposition 
of particle velocities 

fk ^ (/)(V) +£(V), (/)(x,0 = ^^^(x-x ;) 



Here f^m^Xk are velocity, mass and radius-vector of particle number A;; x is coordi- 
nate of the spatial point where velocity is calculated; y/ is localization function; M is 
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the total number of particles in the system. Decomposition (6) can be considered as 
spatial averaging with weight determined by function y/. Note that according to this 
approach thermal component of the velocity /i(x,f) is continuum value (it is deter- 
mined at points between particles) . Several thermal velocities fa (xo , /) , fa (xn , t ) , . . . 
and one continuum velocity ^/\(xo,f) are simultaneously determined in one spa- 
tial point xo. However, formally it does not lead to any contradictions. This type 
of decomposition as well as all mentioned above is not unique. It strongly depends 
on the choice of localization function. In particular, if localization area surrounds 
single atom, then thermal component of the velocity is equal to zero. 

Since the unique decomposition does not exist, the theory should not be based 
on particular method of decomposition. In addition, the results of the theory should 
not qualitatively change with the change of the method. Properties of particular 
methods for averaging are not used in the present paper, unless otherwise stated. Let 

us speak about averaged (f\ and thermal / components of physical value /, which 
are connected by formula (1). 

The second important statement used in the present paper is the long wave ap- 
proximation [2]. The idea of the approximation is the following: an average com- 
ponent of any physical value is assumed to be slowly changing in space on the 
distances of an order of the interatomic distance. Then the average component can 
be considered as a continuum function of a space variable and can be expanded 
into power series with respect to interatomic distance. The resulting series should 
converge rapidly. Exactly this assumption allows to make transition from a discrete 
system to an equivalent continuum. 



3 Kinematics 

Let us use material description of the equivalent continuum. Two configurations of 
continuum and discrete system are considered: reference and actual. For the sake of 
simplicity let us take undeformed configuration of the crystal lattice as the reference 
one. Radius-vectors of equivalent continuum in the reference and actual configura- 
tions are denoted as r and R respectively. Two ways of particles' identification are 
used. On the one hand, the position of the particle is determined by its radius-vector. 
On the other hand, let us use local numbering [7]. Starting with one reference par- 
ticle let us mark all its neighbors by index a. Let us denote a vector connecting the 
reference particle with its neighbor number a as a a . By the definition vectors a a 
have the following property 

&a = — 3-a- (7) 

The same vectors in an actual configuration are represented as a sum of averaged 
component A a and thermal component A a . They can be expressed in terms of vec- 
tors a a and displacements of particles as 

A a =a a +u a -u, A„=u a u. 
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Here u, u a and u, u a are average and thermal components of displacements respec- 
tively. 

The introduced way of particle identification has several properties. Let us use 
the following definition. A physical value determined by the state of one particle is 
called a single-particle value. For example, particle's mass, radius-vector, velocity, 
displacement etc. Let /(rn) be one -particle value that corresponds to the reference 
particle with radius-vector ro in the reference configuration. Denote value /, which 
corresponds to particle number a, as /«(ro). Then the following two designations 
are equivalent 

/«(r ) = /(ro + a«). (8) 

In the framework of this approach the magnitude of physical value / at the point ro 
can be represented in the following equivalent forms 

/(ro) =/a(r -a tt ) =/_ a (r + a a ). (9) 

One can show that for multiplication of two one-particle values / and g the follow- 
ing identities are satisfied 

\fa8](To) = \fg- a ](To + Ra) \f-ag](To) = \fga](To-*a)- (10) 

Hereinafter square brackets mean that all values in them are calculated at the same 
point. Let us also consider values that depend on the state the reference particle and 
its neighbor number a, notably vector connecting two particles and force acting 
between particles. These values have the following property 

h a (r ) = -h- a (ro + a a ). (11) 

If h is a force acting between particles, then equation (11) is a specific form of 
Newton's third law. The following identities are satisfied for one-particle value / 
and the value h 

\f a h a }(r ) = -[fh- a }(r + a a ) \fh- a ](r ) = -\f a h a ](r -a a ). (12) 

Let us consider kinematics of the discrete system in the long wave approxima- 
tion. It is assumed that average values of particle radius-vectors are identical to 
positions of corresponding points of continuum media. Thus if some particle has 
radius-vector r in the reference configuration, then the average value of the radius- 
vector in actual configuration is equal to R(r). The average position of its neighbor 
number a is determined by vector R(r + a«). Then one can show that vectors A a 
and a a , connecting the particles, are related by the following formula 

A„=R(r|a„)-R(r)«a„-VR, (13) 

o 

where V is nabla-operator in the reference configuration. Here the long wave ap- 
proximation was used, which allows to leave the first order terms only. One can 
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see that expression (13) is similar to the formula used in continuum mechanics that 
connects vectors dr and dR. Using equation (13) one can derive relations between 
vectors A a , a a and measures of deformation used in nonlinear theory of elastic- 
ity [14]. For example, the following identity fulfills for Cauchy-Green measure G 

A a = a a a a --G, G = (V R) • (R V). (14) 



4 Equation of momentum balance 

Let us obtain the equation of motion for the equivalent continuum. Thereto let us 
write down the equation of motion for the reference particle and use the decompo- 
sition of motions 

mu = WF«(A a +A a )Y mu = £F a (A a +A a ), (15) 

a a 

where F a is a force acting on the reference particle from its neighbor a; m is par- 
ticle's mass. The first equation from (15) describes slow motion of the system. The 
motion can be considered as motion of continuum media. The second equation de- 
scribes thermal oscillations. One can see that both equations are coupled via the 
argument of the force F a . However, if the dependence of the force on the distance 
between particles is linear, then equations become independent. It reflects the well- 
known fact that harmonic models can not describe coupled thermo-mechanical ef- 
fects such as thermal expansion [13]. Let us conduct the following transformations 
in the first equation from (15). 

OT (i = £(F a ) = £(f_ b ) = i£(F a +F_A (16) 

Force F a satisfies Newton's third law, i.e. F a (r — a«) = — F_ a (r) (see formula (1 1)). 
Averaging this expression and using long wave approximation one obtains 

(F_„)(r) w -^F a )(r)+a„- V (F a )(r) (17) 

Substituting formula (17) into equation (16) and dividing both parts by volume of 
elementary cell in the reference configuration Vo one obtains 

|«=v.(^E.„(f„)) 

Let us compare formula (18) with equation of motion for continuum in Piola's 
form [14]. 

Poii=V-P. (19) 
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where P is Piola stress tensor, po is a density in the reference configuration. Com- 
paring equations (18) and (19) one can deduce that 

P = ^X>( F °}> po = ^ (20) 

Strictly speaking the first formula from (20) is satisfied only with accuracy of tensor 
with zero divergency. This tensor corresponds to some equilibrium stress field in the 
crystal. 

Let us conduct the same derivations in an actual configuration. Equation of mo- 
tion for the particle has form (16). Let us rewrite formula (17) in an actual configu- 
ration. 

F«(R-A a ) = — F_ a (R) ^F_ a )(R)«-^F a )(R)+A a -V^F a )(R). 

(21) 

Here it was used that in the long wave approximation A_ a ss — A a . Substituting ex- 
pression (21) into equation (16) and dividing both parts by the volume of elementary 
cell in the actual configuration V one obtains 

5,_^Ei..v(F.>. (22) 

Let us conduct the following transformations in the right side of formula (22) 

^EA a .v(F a }=V. (^EA a (F a )) - J> (±Aa) (F«) (23) 

The second term in the right side of the given equation can be written down in the 
following form using equation (13) 



RV •a a (F a )=0, (24) 



where Piola's identity V • (iL = was used (see, for example, [14]). Then 
equation of motion (22) has the following form 

The requirement of equivalence of discrete and continuum systems leads to the fol- 
lowing expressions for Cauchy stress tensor and density in the actual configuration 

1 / \ Til 
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If thermal motion is not taken into account, then expression (26) coincides with 
expressions derived in papers [7, 20]. 

It is known that Cauchy stress tensor is symmetrical in systems without moment 
interactions. Let us consider tensor T determined by formula (26). Force F a can be 
represented as 

F a = -4> a ((A a +A a ) 2 )(A a +A a ), 0(A 2 ) d =l f -^l (27) 
Substituting the given expression into formula (26) one obtains 

T = -^L( «) A « A «-^E A «(^ A «)- ( 28 > 

The first tensor in the right side of formula (28) is symmetrical indeed. However 
the symmetry of the second tensor is not evident. Further it will be shown that 
antisymmetrical part of this tensor is small with respect to symmetrical part. 



5 Equation of angular momentum balance 

It is known from continuum mechanics [15] that the symmetry of Cauchy stress 
tensor follows from equation of angular momentum balance for elementary volume. 
In the discrete case elementary cell plays the role of elementary volume. Let us 
write down the averaged equation of angular momentum balance for elementary 
cell (moments are calculated with respect to the center of the cell determined by 
vector R). 

m/u x u)'= (u x £F„) = -£( A « x F„) + ^(s a x F a ). (29) 
a a a 

Transforming the second term in the right side of the given equation using the long 
wave approximation one obtains 

(u a x F«)(R) = -(u x F_ a ^(R + A«) « -(u x F_ a ) -A a • V^5 x F_ a ). 

(30) 

Let us substitute the result into equation (29) and resolve it with respect to ^A a x F a 
i^/A a xF a \ = ^£ A « x /<£aA a \ = -i^A„ • v/SxF_„\ -m/uxuV 

Z a I L a \ / Z a 

(31) 

Using expression (28) for the stress tensor one can transform formula (31) to the 
following form 
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E-xT A = ^ 7 ^A ce -V^5xF_ a )+p^Sxiiy (32) 

Here A denotes antisymmetrical part of the tensor. One can see that if there is no 
thermal motion, then T A = 0. Let us show that in general case T A is small in compar- 
ison with T . The first term in formula (32) is small in the long wave approximation. 
Let us show that if operator (\ contains spatial averaging, then the second term is 
also small. Consider the following identity 

p(uxu\ = p/(R+u)xuY. (33) 

The right side of formula (33) is a derivative of angular momentum, which cor- 
responds to thermal motion. Angular momentum is calculated with respect to the 
origin of coordinates. Let the averaging operator include spatial averaging over 
significantly large volume and let us assume that thermal motion does not lead to 
macroscopic rotation of the volume. Then expressions (33) are equal to zero. As a 
result z A has the same order as terms that were neglected in long wave assumption. 
Consequently, tensor (26) can be considered as approximately symmetrical. 

Thus averaging operator proposed above can not be arbitrary. It should include 
spatial averaging. Otherwise one can not prove the symmetry of tensor T determined 
by formula (26). 



6 Equation of energy balance 

For the sake of simplicity let us assume that volumetrical forces and volumetrical 
heat sources are equal to zero. Derivations are carried out in the reference configu- 
ration. In this case averaged specific total energy per volume Vq has the following 
form 

Po^ = \po( (u + u) 2 ) + ^ £(n(A„ + A«) ) , (34) 

where S is particle's total energy divided by the mass, i.e. discrete analog for mass 
density of the energy. Let us introduce the following designations 

I i /. 2 \ 1 _ / \ (35) 

pnJr=-p u ' Po^ = 2po(u ^ + _£^n(A a +A a )^. 

Values and °i/ correspond to mass densities of macroscopic kinetic and inter- 
nal energies. Calculating derivatives of kinetic and potential energies taking into 
account formulas (15), (20) one can obtain 
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l 2 P °7t (" 2+ <") 2 ) = Pb (*■*+(*■£)) = (V-P) -u+^i:(u.F a ), 

(36) 

_-E(u(A a+ i a ))--i;(F a .(A a+ A a ): 

where formula F a = — — was used. Let us conduct the following transformations 
dA a 

L(f«) ■ A„ = E( F «) ' (u(r + a«) -ii(r)) w £aa(F a ) ■ li V (37) 
o a a 

Summarizing expressions (36) and taking into account formulas (20), (37) one ob- 
tains the expression for derivative of the total energy with respect to time 

Po ,#=V-(P-u) + ^-£(eV (u a +u)). (38) 

Using equation of momentum balance one can show that po Jf? = -P^ • u. Sub- 
stituting this expression into equation (38) one obtains 

Po# = P--(uv) + 2^£(F«-(u ce +u)y (39) 

Comparing the last expression with energy balance equation for a continuum me- 
dia [15] one can conclude that expression for divergency of heat flux in the reference 
configuration h has form 

Let us represent the right side of this expression in the form of divergency. Using 
the first identity form (12) in the right side of formula (40) one obtains 



F a ■ u a ^ (r) = - (F_„ • u )(r + a a ), (F- a ■ u) (r) = - ( F a • u a j (r - a„). 

(41) 

Substituting formulas (41) into formula (40) and applying the long wave approxi- 
mation one can obtain 

V * =V • (-^Ea«(F a -3 a )) =V • ■«)) • (42) 

Using this expression one can write down three representations for heat flux in the 
reference configuration 1 



Note that heat flux is determined with the accuracy of vector with zero divergency. 
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(43) 

Analogous formulas can be obtained for heat flux in the actual configuration H 

H = -^E A «(F« • (U«+U)) = -^E A «(F« -Ua) = -^E A «(F« 

(44) 

Here formula H = y (iL V^j ■ h relating heat fluxes in different configurations was 
used [6]. Note that different expressions for h and H in formulas (43), (44) are equal 
with accuracy of terms, which were neglected in the long wave approximation. 



7 Constitutive relations for stress tensor and heat flux 

Expressions (26), (35) connecting micro- and macro-parameters allow to derive 
nonlinear constitutive relations (equations of state) for thermo-elastic behavior of 
the crystal. This problem is considered in detail in works [7, 9]. Only main ideas 
and results are shown below. The relation between "cold" component of Cauchy 
stress 2 To and Cauchy-Green measure of deformation can be obtained substituting 
formulas (13), (27) into formula (26) (see [7] for details). 

t = ^=(RV)- (£<P(a a a --G)a a a J -(VR), (45) 

2Voy|G| \a / 

where |G| is a determinant of tensor G. Equations of state connecting thermal com- 

def 

ponent of Cauchy stress Zj = T — To with thermal energy were obtained in paper [9] . 
The expansion of Cauchy stress and internal energy with respect to A a was con- 
ducted. In particular, in the first approximation the following system was obtained 

Tt = -^tL [24>«A a EA + 4>' ce A ce A ce E + 2^A ce A ce A ce A ce ] • • (a«A« 

(46) 

Ut = ~I [* a E + 24>' a A a A a ] • • (a A ), ' 



2 i—i I " " T o"""nj y a a / ^ a — d(A 2 )" 

Here Ut is a thermal energy per unit volume. It was assumed that ^A a A a 

2 K 2 E, K 2 = (j^a^ m or der to close system (46). In the framework of the as- 
sumption system (46) takes form 

1 r g M(rf + 2)g B + 2^ B )A B A B (4?) 



2 Stress in the crystal in the absence of thermal motion. 
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Expression (47) is generalized Mie-Gruneisen equation, where F is tensor Gruneisen 
coefficient. Note that more accurate equations of state can be obtained leaving 
higher order terms in expansions (46). For further details about the approach for 
derivation of equations of state see paper [9] . 

Let us consider propagation of small thermal disturbances. Assume that the am- 
plitude of thermal oscillations is small in comparison with interatomic distance. 
Expanding expression for heat flux (44) with respect to A a and leaving terms of 
order of A^ only one obtains 



^£(0(A 2 JA a E + 2<£'(A 2 JA a A C( A C( ) ■■(A a fi a ), ^ = 



Expression (48) is satisfied for arbitrary nonlinear elastic deformations. Let us con- 
sider the case when discrete system is free from internal mechanical loads and 
constraints. In this case deformations are caused by thermal expansion only. Then 
the reference and actual configurations approximately coincide. Linearizing expres- 
sion (48) assuming that A a rj a a one obtains 

; 3 c a ..(A a s a )=x 3 

a a 



H = £ 3 C a • • ( A a u a ) = £ 3 C« • • ( A„5) , (49) 



where 3 C a = {®{a\i)&aR + {a 2 a )& a & a &a). Let us represent the expression 



2V 

for heat flux in the form of divergency 



H = £ 3 C a ■■(A a u a ) w V- (^£a a 3 C a ■ -(SuV) -£ 3 C a 



S 

UU tt ) . (50) 



Here the following identity was used La 3 C a = 0. From formula (50) it follows that, 
in contrast to classical Fourier law, the heat flux depends on the set of symmetrical 

tensors {^uu^), ^uu a ^ . Let us try to connect heat flux with temperature. Classical 
ideal gas definition of temperature is used 

dkT = tn(Z 2 y (51) 

where k is Boltsman constant. Equation (51) can be transformed taking into account 
equation of motion (15) 

dkT = m(u • uY - (£F« • u) « m(u ■ u\+ £ (<PE + 2<P'a a a a ) ■ ■ (uA a ). (52) 



Here the expansion into series with respect to A a is carried out. The second order 
terms are leaved only. Using the definition of tensor 3 C« let us write down the 
resulting system for connection between heat flux and temperature. 
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According to system (53) the thermal state at the given point is determined by sym- 
metrical tensors 3 ^uu^), (uua^ . In general, these tensors are independent. There- 
fore system (53) is not closed. 



8 Concluding remarks 

The generalization of the approach proposed in [7] that allows to carry out the trans- 
formation from discrete system to equivalent continuum was presented. Two main 
principles were used for the transformation: the decomposition of particles' motions 
into continuum and thermal parts, and the long wave assumption [2]. The review 
of different methods for decomposition was given. It was shown that all of them 
contain uncertain parameters. Therefore, the result of decomposition is principally 
nonunique. Thus, one can conclude that derivations should not be based on any 
specific decomposition type. The connection between kinematics of discrete sys- 
tem and kinematics equivalent continuum was analyzed. Equivalent Cauchy-Green 
measure of deformation for discrete system was introduced. The transition form 
single particle's equation of motion to equation of motion for equivalent continuum 
was carried out. Expressions connecting Cauchy and Piola stress tensors with pa- 
rameters of the discrete system were derived. It was shown that discrete analog of 
Cauchy stress tensor can be non-symmetrical. Spatial averaging is necessary for the 
symmetry of this tensor. Thus, averaging operator cannot be arbitrary and should 
contains spatial averaging. The energy balance equation for discrete system was 
considered. The equation was transformed to the form similar to energy balance 
equation for a continuum system. As a result, the expression connecting heat flux 
with parameters of discrete system was obtained. Propagation of small thermal dis- 
turbances in undeformed crystal was analyzed. It was shown that thermal state at the 

point is determined by the set of independent symmetrical tensors ^uu^, (uua^ ■ 
This fact does not allow to connect heat flux with temperature. Equation of state in 
generalized Mie-Gruneisen form connecting Cauchy stress tensor with deformation 
gradient and thermal energy is obtained from microscopic considerations. 
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3 Tensors ^uu„j can be represented via time derivatives of tensors ^uu a ^) using long wave 
approximation. 
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